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Microscopic calculations based on realistic nuclear hamiltonians, while yielding accurate results 
for the energies of the ground and low-lying excited states of nuclei with A < 12, fail to reproduce 
the empirical equilibrium properties of nuclear matter, that are known to be significantly affected 
by three-nucleon forces. We discuss a scheme suitable to construct a density-dependent two-nucleon 
potential, in which the effects of n-particle interactions can be included by integrating out the degrees 
of freedom of (n — 2)-nucleons. Our approach, based on the formalism of correlated basis function 
and state-of-the-art models of the two- and three-nucleon potentials, leads to an effective interaction 
that can be easily employed in nuclear matter calculations, yielding results in good agreement with 
those obtained from the underlying three-body potential. 



PACS numbers: 21.30.Fe, 21.45.Ff, 21.65.-f 

I. INTRODUCTION 

The results of ab initio microscopic calculations consis- 
tently suggest that realistic nuclear hamiltonians, includ- 
ing both two- and three-nucleon potentials, while provid- 
ing a quantitative account of the energies of the ground 
and low- lying excited states of nuclei with A < 12 [HQ, 
fail to explain the empirical equilibrium properties of nu- 
clear matter. This problem can be largely ascribed to 
the uncertainties associated with the description of three- 
nucleon interactions, whose contribution turns out to be 
significant. 

A signal of the limitations of the commonly employed 
three-nucleon potential models (e.g. the Urbana IX 
model of Ref. has been recently provided by the 

authors of Ref. [J], who carried out a study of sym- 
metric nuclear matter within the Auxiliary Field Diffu- 
sion Monte Carlo (AFDMC) approach. Their results, ob- 
tained using a truncated version Hof the state-of-the-art 
nucleon-nucleon potential of Ref. [fj , show that AFDMC 
simulations do not lead to an increase of the binding 
energy predicted by Fermi-Hyper-Nettcd-Chain (FHNC) 
and Brueckner-Hartree-Fock (BHF) calculations 0]. 

Different three-nucleon potential models 0, H| , yielding 
similar results when applied to the calculation of nuclear 
properties, predict sizably different equations of state 
(EoS) of pure neutron matter at zero temperature and 
densities exceeding the nuclear matter saturation den- 
sity, po = 0.16 fm -3 [9(. In this region, the three-nucleon 
force contribution to the binding energy becomes very 
large, the ratio between the potential energies associated 
with two- and the three-body interactions being ~ 20% 
at density p ~ 2po (see, e.g. Ref. [Io|). The size of 



the three-body force contribution suggests that, at large 
p, interactions involving four or more nuclcons may also 
play an important role, and should be taken into account. 

In view of the severe difficulties involved in the imple- 
mentation of the existing models of three-nucleon inter- 
actions in many-body calculations, the explicit inclusion 
of four- and more-body potentials does not appear to 
be a viable option. In this paper we follow a different 
strategy, somewhat along the line of the Three- Nuclcon- 
Interaction (TNI) model proposed by Lagaris and Pand- 
haripande [ll| and Friedman and Pandharipande [l2T | in 
the 1980s. 

The authors of Refs. [ill, suggested that the main 
effects of three- and many-nucleon forces can be taken 
into account through an effective, density-dependent two- 
nucleon potential. However, they adopted a purely phe- 
nomenological procedure, lacking a clearcut interpreta- 
tion based on the the analysis of many-nucleon interac- 
tions at microscopic level. 

The TNI potential consists of two density-dependent 
functions involving three free parameters, whose values 
were determined through a fit of the saturation density, 
binding energy per nucleon and compressibility of sym- 
metric nuclear matter (SNM) , obtained from FHNC vari- 
ational calculations. The numerical values of the three 
model parameters resulting from recent calculations per- 
formed by using AFDMC simulations turn out to be only 
marginally different from those of the original TNI po- 
tential [l]|. 

The TNI potential has been successfully applied to ob- 
tain a variety of nuclear matter properties, such as the 
nucleon momentum distribution [14j | , the linear response 
[IH [Hj|, and the Green's function 0, El], within the 
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Correlated Basis Function (CBF) approach (for a review 
of CBF theory and its applications, see Ref.[19| and ref- 
erences therein). 

The strategy based on the development of two-body 
density-dependent potentials has been later abandoned, 
because their application to the study of finite nuclei 
involves a great deal of complication, mainly stemming 
from the breakdown of translation invariance. While in 
uniform matter the density is constant and the expansion 
of the effective potential in powers of p is straightforward, 
in nuclei different powers of the density correspond to dif- 
ferent operators, whose treatment is highly non trivial. 

However, the recent developments in numerical meth- 
ods for light nuclei seem to indicate that the above diffi- 
culties may turn out to be much less severe then those im- 
plied in the modeling of explicit many-body forces and, 
even more, in their use in ab initio nuclear calculations. 

In view of the observation, based on a variety of ex- 
perimental evidence |2(| HH , that short range nucleon- 
nucleon (NN) correlations are a fundamental feature of 
nuclear structure, the description of nuclear dynamics in 
terms of interactions derived in coordinate space, like the 
Urbana-Argonne models, appears to be the most appro- 
priate, for both conceptual and technical reasons. 

First of all, correlations between nucleons are predomi- 
nantly of spatial nature, in analogy with what happens in 
all known strongly correlated systems. In addition, one 
needs to clearly distinguish the effects due to the short- 
range repulsion from those due to relativity Finally, 
quantum Monte Carlo methods have serious difficulties 
in dealing with highly non local interactions. For all the 
above reasons we stick to two-body density-dependent 
potentials of the Urbana-Argonne type. 

Our approach is based on the tenet that n-body poten- 
tials (n > 3) can be replaced by an effective two-nucleon 
potential, obtained through an average over the degrees 
of freedom of n — 2 particles. Hence, the effective po- 
tential can be written as a sum of contributions ordered 
according to powers of density, the p-th order term being 
associated with (p + 2)-nucleon forces. 

Obviously, such an approach requires that the average 
be carried out using a formalism suitable to account for 
the full complexity of nuclear dynamics. Our results show 
that, in doing such reduction, of great importance is the 
proper inclusion of both dynamical and statistical NN 
correlations, whose effects on many nuclear observablcs 
have been found to be large [2(| Hl( . 

In this work, we use CBF and the Fantoni-Rosati (FR) 
cluster expansion formalism [l9|, [22|, HU to perform the 
calculation of the terms linear in density of the effective 
potential, arising from the irreducible three-nucleon in- 
teractions modeled by the UIX potential. 

It should be noticed that our approach significantly 
improves on the TNI model, as the resulting potential is 
obtained from a realistic microscopic three-nucleon force, 
which provides an accurate description of the properties 
of light nuclei. 

While being the first step on a long road, the results 



discussed in this paper are valuable in their own right, 
as the effective potential can be easily implemented in 
the AFDMC computational scheme to obtain the EoS of 
SNM. Similar calculations using the UIX potential are 
not yet possible, due to the complexities arising from the 
commutator term. In addition, the density-dependent 
potential can be used to include the effects of three- 
nucleon interactions in the calculation of the nucleon- 
nucleon scattering cross section in the nuclear medium. 
The knowledge of this quantity is required to obtain a 
number of nuclear matter properties of astrophysical in- 
terest, ranging from the transport coefficients to the neu- 
trino emission rates [13, [H[ . 

In Section [H] we discuss the main features of the exist- 
ing theoretical models of the three-nucleon force, while 
Section IIIII is devoted to a brief review of the many- 
body approach based on CBF and the cluster expansion 
technique. In Section IIVI we describe the derivation of 
the density-dependent interaction, pointing out the role 
of dynamical and statistical correlation effects. In Sec- 
tion [V] we compare the energy per particle of nuclear 
matter obtained from the effective potential to that re- 
sulting from highly refined calculations, carried out using 
the Argonne v' 6 [26j j and v' s (H nucleon-nucleon potentials 
and the Urbana IX three-nucleon potential p|. Finally, 
in Section I VII we summarize our findings and state the 
conclusions. 



II. THREE NUCLEON FORCES 

Nuclear many-body theory (NMBT) is based on the 
assumption that nuclei can be described in terms of point 
like nucleons of mass m, whose dynamics are dictated by 
the hamiltonian 

i j>i k>j>i 

Before describing the three nucleon potential Vijk, let 
us discuss the vg, two-body potential model, that will be 
used throughout the paper. It is given by 

fltf=£^(»-tf)Q&. (2) 
■p=i 

where 

off 1 - 8 = (1,^,%,^, • a tj ) ® (I,-*,-) . (3) 

In the above equation, er^ = er^ • (Tj and = t± ■ Tj, 
where er 2 and Tj are Pauli matrices acting on the spin or 
isospin of the i-th, while 

S tJ = iffafof = (3fgrg. - <P>?of , (4) 

with a, (3 = 1, 2, 3, is the tensor operator, Ly- is the 
relative angular momentum 

Uj = ^(ri-rj) x (V,- Vj) (5) 
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and Sij is the total spin of the pair 

s y = ^( <T » + cr j)- ( 6 ) 

Such potentials have exactly the same form as the first 
eight components of the state-of-the-art Argonne Vig po- 
tential Q . We will be using the so called Argonne v' 8 and 
Argonne Wg potentials, which are not simple truncations 
of the Argonne v\g potential, but rather rcprojections 

The Argonne v' s potential is obtained by refitting the 
scattering data in such a way that all S and P partial 
waves as well as the 3 Z?i wave and its coupling to 3 Si 
are reproduced equally well as in Argonne vi&. In all 
light nuclei and nuclear matter calculations the results 
obtained with the v' s are very close to those obtained 
with the full vis, and the difference v\g — v' s can be safely 
treated perturbatively. 

The Argonne v' 6 is not just a truncation of v' s , as 
the radial functions associated with the first six oper- 
ators are adjusted to preserve the deuteron binding en- 
ergy. Our interest in this potential is mostly due to the 
fact that AFDMC simulations of nuclei and nuclear mat- 
ter can be performed most accurately with ug-type of 
two-body interactions. Work to include the spin-orbit 
terms in AFDMC calculations is in progress. On the 
other hand we need to check the accuracy of our pro- 
posed density-dependent reduction with both FHNC and 
AFDMC many-body methods before proceeding to the 
construction of a realistic two-body density-dependent 
model potential and comparing with experimental data. 

It is well known that using a nuclear Hamiltonian in- 
cluding only two-nuclcon interactions leads to the under- 
binding of light nuclei and overestimating the equilibrium 
density of nuclear matter. Hence, the contribution of 
three-nucleon interactions must necessarily be taken into 
account, by adding to the Hamiltonian the corresponding 
potential, e.g. the widely used Urbana IX (UIX) Q. 

The potential of Ref. Q consists of two terms. The 
attractive two-pion (2w) exchange interaction V 2n turns 
out to be helpful in fixing the problem of light nuclei, but 
makes the nuclear matter energy worse. The purely phc- 
nomenological repulsive term V R prevents nuclear matter 
from being overbound at large density. 

The V 2lT term was first introduced by Fujita and 
Miyazawa [27j to describe the process whereby two pi- 
ons are exchanged among nuclcons and a A resonance is 
excited in the intermediate state, as shown in the Feyn- 
man diagram of Fig. [T] It can be conveniently written in 
the form 

v 2n = E = E ( ? ) 

cyclic cyclic 

where 

0?23 =^({A 12 ,A 23 }{t 12 ,t 23 } 

+ ^[A i2 , A > 23 ][ti 2 ,t 23 ]^ (8) 
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Figure 1. Feynman Diagram associated with the Fujita 
Miyazawa three-nucleon potential term. 



and 

Xij = Y(m„r)<Tij + T(mnr)Sij . (9) 

The radial functions associated with the spin and tensor 
components read 

Y(x) = —d Y (x) (10) 

x 

T{x)=(l + - + \)Y{x)^{x), (11) 

V X X* J 

while the £(x) are short-range cutoff functions defined by 
t Y (x)=t T {x) = l-e- cx \ (12) 

In the UIX model, the cutoff parameter is kept fixed at 
c = 2.1 fm -2 , the same value as in the cutoff functions 
appearing in the one-pion exchange term of the Argonne 
V\s two-body potential. On the other hand, A 27r is varied 
to fit the observed binding energies of 3 H and 4 He. The 
three-nucleon interaction depends on the choice of the 
NN potential; for example, using the Argonne v\% model 
one gets A 27r = -0.0293 MeV. 

The repulsive term V R is spin-isospin independent and 
can be written in the simple form 

yR = E V £ = f«E T 2 {m^ t] )T 2 {m^ ]k ) , (13) 

cyclic cyclic 

with T(x) defined in Eq. (fTTj) . The strength Uq, adjusted 
to reproduce the empirical nuclear matter saturation den- 
sity, is U = 0.0048 MeV with v w . 

The two parameters A 27r and Uq have different values 
for v' 8 and v' e . We disregard such differences in the 
present analysis, mostly aimed at testing the quality of 
our density-dependent reduction of the UIX three-body 
potential, rather than reproducing empirical data. 
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III. FORMALISM 

A. Correlated basis theories 

One of the most prominent features of the nucleon- 
nucleon (NN) interaction is the presence of a repulsive 
core, giving rise to strong correlations that cannot be 
taken into account within the independent particle pic- 
ture. 

This problem has long been recognized, and was clearly 
pointed out by Blatt and Wcisskopf over fifty years ago. 
In their classic Nuclear Physics book, first published in 
1952, they warn the reader that "the limitation of any 
independent particle model lies in its inability to encom- 
pass the correlation between the positions and spins of 
the various particles in the system" (29|. 

Let us consider uniform nuclear matter, defined as 
a translationally invariant system of protons and neu- 
trons, in which the electromagnetic interaction is turned 
off. In the absence of interactions, such a system can 
be described as a Fermi gas at zero temperature, and 
its ground state wave function reduces to the antisym- 
metrized product (Slater determinant) of orbitals associ- 
ated with the single particle states belonging to the Fermi 
sea: 



$(xi 



,x A ) = A[4>„^{xx) ...<j) nA {x A )'. 



(14) 



with 



and 



= ^ki^T^r*) = ^(r^Xo-i^Ti , (15) 



(16) 



In the above equations, i7 is the normalization volume, 
Xa-i and rj Ti are Pauli spinors, describing the nucleon spin 
and isospin and |ki| < kp = (6n 2 p/v) 1 ^ 3 . Here kp is 
the Fermi momentum, while p and v denote the density 
and the degeneracy of the momentum eigenstates [v =2, 
4 pure neutron matter (PNM) and symmetric nuclear 
matter (SNM), respectively]. The generalized coordinate 
Xi = {fi, c;, Ti \ represents both the position and the spin- 
isospin variables of the i-th nucleon, while n.i denotes 
the set of quantum numbers specifying the single particle 
state. 

The antisymmetrization operator A can be written in 
the form 

A = l-J^P ij + ( p *i p H + p ikPkj) + ■ ■ ■ , (17) 

i<j i<j<k 

where 

Pij = t(1 + <Tij)(l + Ty) expH(ki - kj) • ry] (18) 

is the two-particle exchange operator, defined by the re- 
lation 



Note that, as shown by Eq. (|T8|) . the exchange operators 
act on both the radial and spin-isospin components of 
the nucleon wave function. 

Due to the strong repulsive core, the matrix elements 
of Vij between eigenstates of the non interacting system 



^k 1 ,cr 1 ,T 1 ,</ , k 2 ,(T 2 ,T 2 , l^l'/'kicriTi^ka 



(20) 



turn out to be very large, or even divergent if the core of 
the NN potential is infinite. As a consequence, perturba- 
tive calculations carried out using the bare NN potential 
and the Fermi gas basis states are unavoidably plagued 
by lack of convergence. 

To circumvent this problem, one can follow two differ- 
ent strategies, leading to either G-matrix or CBF per- 
turbation theory. Within the former approach, the bare 
potential is replaced by a well behaved effective in- 
teractions, the so called G-matrix, which is obtained by 
summing up the series of particle-particle ladder dia- 
grams. In the second approach, nonperturbative effects 
are handled through a change of basis functions. 

Correlated basis theories of Fermi liquids [ill l22l l33l . l34j 
are a natural extension of variational approaches in which 
the trial ground state wave function is written in the form 



l*o> 



(^oli^N) 1 / 2 



(21) 



In the above equation, F is a correlation operator, whose 
structure reflects the complexity of the nucleon-nucleon 
potential (3l| : 



with 



F = S J] Fij , 

J>»=1 



p=l 



(22) 



(23) 



Note that the symmctrization operator S is needed to 
fulfill the requirement of antisymmetrization of the state 

basis (CB) is defined as 



l^n), since, in general, [Of- , Of fe ] ^ 0. The correlated 



l*n> - 



F\<S>n) 



;$„|FtF|$„)V2 



(24) 



Ptj4>ni{Xi)<t>nAXj) = 4>n,(Xj)4>nAXi) 



(19) 



where is a generic n particle - n hole Fermi gas state. 
The CB states are normalized but not orthogonal to each 
other. They have been used within non orthogonal per- 
turbation theory j33l [35| to study various properties of 
quantum liquids. An exhaustive analysis of the conver- 
gence properties of this perturbation scheme has never 
been carried out, but the truncation of the series at a 
given perturbative order is known to lead to nonorthogo- 
nality spuriositics, whose effects are not always negligible. 
A much safer and efficient procedure, in which one first 
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orthogonalizes the CB states by using a combination of 
Schmidt and Lowdin transformations and then uses stan- 
dard perturbation theory, has been proposed by Fantoni 
and Pandharipande [34[. 

The radial functions f p (rij), appearing in the defini- 
tion of the correlation operator are determined by the 
minimization of the energy expectation value 

E v = (*o|#|*o>, (25) 

which provides an upper bound to the true ground state 
energy Eq. In principle, that can be done by solving 
the Euler equations resulting from the functional mini- 
mization of Ey with respect to the correlation functions 
f p (rij), in analogy with what has been done in Jastrow 
theory of liquid 3 He. However, the presence of the spin- 
isospin dependent correlation operators and their non- 
commutativity makes the application of this procedure 
to nucleonic systems almost prohibitive. In this case 
the functional minimization can be carried in a more 
straightforward fashion on the lowest order cluster con- 
tribution to Ey, paying the price of introducing proper 
constraints and the associated variational parameters, as 
discussed below. 



B. Cluster expansion 

In CBF theories the calculation of Ev is carried out 
by i) expanding the r.h.s. of cq. ([25)1 in powers of proper 
expanding functions that vanish in uncorrelated matter 
and ii) summing up the main series of the resulting clus- 
ter terms by solving a set of coupled integral equations. 
The FR cluster expansion [23| has been derived to accom- 
plish the first of these two steps for the case of Jastrow 
correlated models. It has been obtained through a gener- 
alization of the concepts underlying the Mayer expansion 
scheme, originally developed to describe classical liquids 
[30j . to the case of quantum Bose and Fermi systems. In 
this case the expanding quantity is given by 

Kr ij ) = f c {r ij f-l, (26) 

where f c is the only correlation of the Jastrow model. 
Notice that in the calculation of Ev the kinetic energy 
operators V^ 2 also act on the correlation functions, giv- 
ing rise to additional expanding quantities. For the sake 
of simplicity, and since here we are interested in cal- 
culating the expectation values of two- and three-body 
potentials, we will not address this issue. It has been 
proved [23[ that the FR cluster expansion is linked, and 
therefore does not suffer the appearance of infinities in 
the thermodynamic limit. The FR techniques have been 
subsequently extended and extensively used to deal with 
spin-isospin dependent correlation operators, like those 
of Eq. (p2| (3lL I32I ] . In this case, besides the expanding 
function h(rij) of eq. (f26|) one has to also consider the 
following ones: 

2/ c (r u -)/ p>1 M , / J,>1 (^)/ 9>1 (^). (27) 



The cluster terms are most conveniently represented 
by diagrams fl9l l3lj . The diagrammatic representation 
of the above expanding functions is given by the bonds 
displayed in Fig. [SJ in which hij is represented by a 
dashed line, Zffjffj by a single wavy line, and ff^fi^ 1 
by a doubly wavy line . 

• • hij 

i j 




Figure 2. Different kinds of correlation bonds. 

The h(rij) is the largest of all the expanding quantities 
and the cluster terms (and similarly the corresponding 
cluster diagrams) involving these functions (or bonds) 
have to be summed up as massively as possible. This 
can be accomplished by solving the FHNC equations of 
Ref. (36j . which already in their basic form sum up all 
diagrams at all orders, with the exception of the so called 
elementary diagrams. 

The cluster diagrams involving operatorial bonds, like 
those representing the functions given in Eq. (f2"T)) . can- 
not be summed up as massively as the scalar diagrams 
of FHNC type. This is due to the additional complex- 
ity associated with the non commutativity of the spin- 
isospin dependent correlation operators. The most pow- 
erful summation scheme which has been derived so far is 
the so called Single Operator Chain of Rcfs. ((3lL l33[). 
generally denoted as FHNC/SOC approximation. 

For the sake of clarity, in the following we summarize 
the main features of the FR cluster expansion and the 
FHNC/SOC approximation, extensively reviewed in [l9|, 

Sl|. 

Let us consider the expectation value of the NN po- 
tential. Exploiting the symmetry properties of the wave 
function it can be written in the form 

($) = (** \FH 12 F\$ ) 
2 (^|FtF|$ ) 

Numerator and denominator of the above equation are 
expanded in powers of the functions defined above. The 
expansion of F^OF, with O = «i2,l for the numera- 
tor and the denominator, leads to series of terms, say 
Xn N ' D \ where the labels N and D stand for numerator 
and denominator, respectively, each characterized by the 
number n of correlated nucleons, i.e. those appearing 
in the argument of the expanding functions. Integrating 
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such terms over the variables of the remaining uncorre- 
cted nucleons amounts to multiplying X n by the n-body 
Fermi distribution operator gp(l, . . . ,n). Consider for in- 
stance one of the cluster terms of the numerator, whose 
structure is given by 



where the Slater function likprij) is given by 



A\ JdX $gX(l,...,n)$ 
/ dX ®*<S> 

dx\ . . . dx n X(l, . . . , n)g F (l, . . . , n) . (29) 



{Xn) 



where dX = dx\ . . . dxA and J dxi stands for integra- 
tion over the coordinate r*j and tracing over the spin and 
isospin variables of the i-th nucleon, and 

. n , A...(A-n + l) 
g F (l,...,n) = 

p n 

JdX<£*<t> { 1 

Note that X(l, . . . , n) can be moved to the left of $g to 
obtain the second line of Eq. (|29|) as, on account of the 
property 

. n , 



dX $*A(l,...,n)$ 



E 



X(l, 



dri . . . dfnl^ixx) . . . <j)* nn {x n )] 



-A[(f> ni (xi) ...(p 



(31) 



one needs to antisymmetrize $o only. 

Summing over the states belonging to the Fermi sea 
for each n$ independently leads to an expression which 
does not depend on the number of particle A. There is no 
violation of the Pauli principle because of the antisym- 
metrization of $o- More specifically, each term of the 
r.h.s. of Eq. ([31]) coming from the antisymmetrization of 
$o is Pauli violating, but the total sum it is not. On the 
other side, the independence of A has very useful conse- 
quences. One of them is that the numerator of (v) can 
be easily recognized as the product of the denominator 
times the sum of linked cluster terms. In addition, the 
FR cluster expansion is exact for any number of particles, 
not just in the thermodynamic limit like, for example, 
the Mayer expansion. This property has been exploited 
in FHNC calculations of finite nuclear systems like nuclei 
[22^ or nucleon confined in periodical box (37| . 

The operatorial n-body Fermi distribution function 
...,n) includes a direct term corresponding to 1 
in Eq. (|17p and a number of exchange terms generated 
according to the algebra of the exchange operators Pij . 

The basic statistical (exchange) correlation is de- 
scribed by the Fermi gas one-body density matrix 



n 

-£(k F rij) XaiVnXijVtj 



(32) 



sin(fci?r,;j) — kpTij coslkpTij) 



(kpnjf 



(33) 



Diagrammatically, the exchange correlation £ij , referred 
to as "exchange bond", is represented by an oriented 
solid line. The -algebra implies that the exchange 
bonds form closed loops which never touch each other. 
If X(l, . . . , n) is made of scalar correlations h{rij) only, 
then all nucleons in a given exchange loop must be in 
the same spin-isospin state and the Fermi distribution 
operators of Eq. (|30|) reduces to the standard Fermi 
gas distribution functions. For example, the two-body 
Fermi distribution function is given by 



1 



gF(rij) = 1 / (k F r i:j ) 



(34) 



As an example, consider the two-body cluster contri- 
bution. From Eq. (|29l) it can be written as 



(X 2 ) = y / ^1^2^(1,2)^(1,2) 



^ [ dndr 2 <f) n M)<t>n 2 (x2yX(l,2) 
x(l-PiaX(n)&, a (z 2 ) (35) 



El 
2 



The sum over the states belonging to Fermi sea implies 
a sum over the spin-isospin quantum numbers, which 
amounts to computing the trace of the spin and isospin 
operators appearing in X(l, 2)(1 — P12). The trace is nor- 
malized to unity, as summation over the momenta leads 
to the appearance of a factor (1/v) in both the direct and 
exchange term. The final result is 

(X*) = \ I d Xl dx 2 X{l, 2)(1- Pi2e 2 (k F r 12 )) . (36) 



1. Diagrammatic rules 

The diagrams consist of dots (vertices) connected by 
different kinds of correlation lines. Open dots represent 
the active (or interacting) particles (1 and 2), while black 
dots are associated with passive particles, i.e. those in the 
medium. Integration over the coordinates of the passive 
particles leads to the appearance of a factor p. 

Active correlations must be treated differently from 
the passive ones, as the components i>f 2 of the two-body 
potential may be singular, thus leading to divergent in- 
tegrals. In the diagrammatic expansion of (v)/A, the 
quantity ^12^12^12 is represented by a thick solid line, 
denoted as "interaction line" and depicted in Fig. [3] 

Exchange lines form closed loops, oriented clockwise or 
counterclockwise, the simplest of which is the two-body 
loop yielding a contribution —£ 2 (kprij)/v. In addition 
to the extra factors coming from the algebra arising from 
the spin-isospin structure of the corresponding cluster 
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term, an n-vertex loop carries a factor (— )(2v)(— l/v) n , 
where — l/i> is associated with each exchange operator £ij 
and —2v is due to the presence of v spin-isospin species 
of the loop, combined with the existence of 2 different 
orientation and to the minus sign coming from the per- 
mutations. The two-body loop is an exception to this 
rule, because there is only one such loop. Therefore, the 
corresponding factor is — \jv rather than —2/ V. 

The correlation bonds of Figs. PEj) cannot be super- 
imposed to each other. They can only be superimposed 
to exchange bonds. 

The allowed diagrams are all linked, as a result of the 
linked cluster property discussed above. 



3 




1 1, P, q 2 

Figure 4. Example of diagram appearing in the cluster ex- 
pansion of {v}. 



O 

1 



<3 fLv p 12 f q 12 ol 2 o[ 2 o 

2 



p, q 



Figure 3. Graphical representation of an interaction line. 



A typical diagram of the FR cluster expansion is 
sketched in Fig. 2] Its contribution to the potential en- 
ergy expectation value (v) is given by 

p 3 

v),„, n m = 3 — 



P \ f 

e 2 (k F r 13 ) 



f c (r 23 )r(r 23 )f l (r 12 ) V P(r 12 )P(r 12 ) 



x -Tr 123 



4Pi3({O r 23: O l 12 }O p 12 
0[ 2 O p 2 {O q 12l O r 23 }) 



(37) 



where the trace Tri23 is carried out in the spin-isospin 
spaces of particles 1, 2 and 3. The factor SI comes from 
the relation J dr\dr 2 = Jl J dri 2 , due to translation in- 
variance, and 3 is a symmetry factor. The four orderings 
appearing on the r.h.s. of Eq. (j37|) correspond to the 
two possible positions of f r (r 23 ), on cither the left- or 
right-hand side of the operator Fi 2 vi 2 F\ 2 . 



2. FHNC/SOC approximation 

All the linked cluster diagrams or sub-diagrams built 
with scalar passive bonds only, with the only exception 
of the so called elementary diagrams, can be summe d up 
in closed form by solving the FHNC equations [ljl |36| . 
The contributions associated with the elementary dia- 
grams can be formally included in the FHNC equations, 
but there is no exact procedure to sum all of them. They 
can only be taken into account approximatively, by ex- 
plicit calculation of the n-point (n > 4) basic structures 
(FHNC/n approximation). However, it is well known 
that in nuclear matter calculations the FHNC approxi- 
mation provides very accurate results. 



On the other hand, diagrams having one or more pas- 
sive operatorial bonds are calculated at leading order 
only. Such an approximation is justified by the observa- 
tion that operatorial correlations are much weaker than 
the scalar ones. Based on this feature, one would be 
tempted to conclude that the leading order amounts to 
dressing the interaction line with all possible FHNC two- 
body distribution functions. This is not true as, besides 
the short range behavior, the intermediate range behav- 
ior of NN correlations also plays an important role that 
needs to be taken into account. In particular, tensor cor- 
relations, and to some extent also exchange correlations, 
have a much longer range than the central ones. 

In order to handle this problem, summing the class 
of chain diagrams turns out to be to be of great impor- 
tance. To see this, let us consider an extreme example 
of a long range correlation, namely a function h(r) that 
heals to zero as a/r 2 , implying that its Fourier transform 
h(k) behaves as (3/k in the long wavelength limit. Chain 
diagrams of /i-bonds are calculated by means of the con- 
volution integral of the various h(rij) in the chain. In 
Fourier space convolution integrals are given by products 
of h. One can easily verify that, in the long wavelength 
limit, any chain diagram is more singular than h. On 
the contrary, the sum of all the chain diagrams has ex- 
actly the same degree of singularity. Hence, summing 
up the series of chain diagrams takes care of long range 
correlations 0- 

The above issue is taken care of by summing up 
the Single Operator Chains (SOC) in the corresponding 
FHNC/SOC approximation [H 132[. SOC are chain di- 
agrams in which any single passive bond of the chain 
has a single operator of the type f c (rij)f p (rij)0 P j or 
—h(rij)£(kprij) x P^, with p < 6, or FHNC-dressed ver- 
sions of them. Note that if a single bond of the chain is of 
the scalar type then the spin trace of the corresponding 
cluster term vanishes, as the Pauli matrices are traceless. 



1 This feature is critical to the calculation of the long wavelength 
limit of the static structure function and the phonon excitations. 
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Then the SOC is the leading order, and at the same time 
it includes the main features of the long range behavior 
of tensor and exchange correlations. 

The calculation of SOC, as that of FHNC chains, is 
based upon the convolution integral of the functions cor- 
responding to two consecutive bonds. Unlike FHNC 
chains, however, the SOC have operatorial bonds. There- 
fore, the basic algorithm is the convolution of two op- 
eratorial correlations having one common point. Let 
us consider two such correlation operators, say X ik = 
£ p =i6* P (m)Of* and Y kj = E p=1 , p y p (r kj )0^. Their 
convolution gives rise to a correlation operator of the 
same algebraic form Z%j = J2 P =i 6 zP { r v)^ij- 



Zij = p \ dx k X lk Yy 



kj 



(38) 



where the functions ff qr 



ellr 



operatorial the convolution vanishes, i.e. £ f 
The explicit expressions of £ pq k 



j - j. depend on the internal angles 
of the triangle ijk. The above equation includes also the 
convolution of the scalar correlations, which is already 
taken into account by the FHNC chain equations. Hence, 
di r . If one of the bonds is scalar and the second is 

rf>lr _ n 

can be found in Refs. 

mm- 

The ordering of the operators within an SOC is imma- 
terial, because the commutator [Oik,O k j} is linear in a k 
and Tfc, and Pauli matrices are traceless. The only or- 
derings that matter are those of passive bonds connected 
to the interacting points 1 or 2. The reason is that the 
interaction line may have up to four operators. There- 
fore, 1 or 2 may be reached by up to five operators and 
one has to take into account the different orderings. The 
underlying spin algebra is lengthy but straightforward, 
and it is given in Ref. [3lj |. As an example, consider the 
cluster diagram of Fig. [5] and the corresponding cluster 
term of Eq. §7$). The two orderings {O r 2dn O l 12 \O p 12 O q 12 
and O l 12 O p 2 {O q 2 , 0^ 3 } give rise to the same trace, which 
in the case oil=p = q = r = 2= a turns out to be 18. 
The full expression of (v)^ Fig m can be easily extracted 
from the energy term W c (de) displayed in Eq. (7.7) of 
Ref. (HJ , and written in terms of the matrices K lpq and 
L lpq and the vector A m , defined as follows 



6 q = V K lpq 6 1 6 P 



L lpq = ±A q K lpq . 



vr.Mrro 1 :,; 

where the + sign applies if 



4" 



o, 



while the — sign applies if 

^iik{O p ij {0%,0^ k }0\ k )=Q. 



(39) 
(40) 
(41) 



The if-matrices are associated with normal orderings, 
like P jO q jOj k O l ik , whereas the L-matrices apply to al- 
ternate orderings, like CF i -0 T - k O\ i O\ k . 

A second important contribution which is included in 
FHNC/SOC approximation is the leading order of the 
vertex corrections. They sum up the contributions of sets 
of subdiagrams which are joined to the basic diagram- 
matic structure in a single point. Therefore, a vertex 
correction dresses the vertex of all the possible reducible 
subdiagrams joined to it. The FHNC equations for the 
full summations of these one-point diagrams are given 
in Ref. p]|. In the FHNC/SOC approximation they are 
taken into account only at the leading order, i.e. includ- 
ing single operator rings (SOR), which are nothing but 
loops of SOC. Vertex corrections play an important role 
for the fulfillment of the sum rules. 



3. Two-body and three-body distribution functions 



The expectation value (|28p. can be conveniently rewrit- 
ten in the form 



W 
A 



p 



dridr 2 vf 2 gf 2 



where 



p 

9l2 



A(A - 1) Tr 12 fdx 3 ... dx A %F^O p 2 F^o 



J dX $*FtF$ 



(42) 



(43) 



are the operatorial components of the two-body distri- 
bution function. 

The FHNC diagrams are divided in 4 separate classes, 
characterized by the type of bonds ending at the vertices 
associated with particles 1 and 2. The different types 
of vertices are denoted "d" for direct, i.e. involving no 
exchange lines, "e" for exchange, i.e. the vertex of an 
exchange loop, and "c" for cyclic, i.e. the vertex of an 
exchange line. Using this notation we can write, 



9 P =9 P dd + 9 P de + 9 P ed + 9l 



(44) 



The two-body distribution functions satisfy the follow- 
ing sum rules 



dn 2 (g c (ri 2 )-i) = -i, 



dri 2 g a (ri 2 ) = -3 , 



p / dri 2 g r7T (n 2 ) = 9. 



(45) 



Note that the above relations also hold true for the dis- 
tribution functions g p of the Fermi gas model, as well as 
for those obtained retaining the f c correlations only. An- 
other sum rule is given by the expectation value of the 
kinetic energy, which can be written in three equivalent 
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forms, known as Pandharipandc-Bethe (PB), Jackson- 
Feenberg (JF) and Clark- Westhaus (CW). In an exact 
calculation they would all give the same results. Numer- 
ical differences between them gauge the degree of accu- 
racy of the approximations employed in the calculation. 

The generalization of Eq. (|28j) to the case of a three- 
body potential, e.g. the UIX model, reads 



(V) 



A! ($S|Ftf 123 F|$ ) 



(A-3)!3! 



($;|FTF|<I>o) 



(46) 



As for the case of the two-body distribution functions 
g\ 2 , it is useful to define three-body distribution func- 
tions gf 23 , reflecting the operatorial structure of V123 
given in Eqs. © and (fT3"|) . Let us write V123 as a sum 
of spin-isospin three-body operators multiplied by scalar 
functions, depending on the relative distances only 



r, 



12:1 



o p 

23^123 



(47) 





Figure 5. Diagrammatic representation of Eq. (|51[) : the two- 
body density-dependent potential includes the effects of both 
the bare three-body potential and the correlation and ex- 
change lines. While g2 dresses the line joining particles 1 
and 2, the dressing being depicted by a line with a big bubble 
in the middle, gs dresses the lines 1 — 2, 1 — 3, and 2 — 3. 



From Eqs. ([8]) and (fT3|) it follows that the sum of the 
above equation involves 19 operators. The expectation 
value of £123 can be written as 



~T = T\ P 2 Y1 [ dr 12 dr 13^23 5l23. 



(48) 



with 



v 

5l23 



A\ Tri 23 Jdx 4 ... dx A <f>lF^O p 123 F$ 



(A -3)1 



p 3 j dX $jFtF$ 



(49) 

In Ref . [38( , the above expectation value has been com- 
puted in FHNC/SOC. The cluster expansion and the cor- 
responding diagrammatic rules of the cluster diagrams 
are very similar to those outlined in the case of the two- 
body potential, with the only difference of three external 
points, instead of two, all linked by interaction lines. 



IV. DERIVATION OF THE EFFECTIVE 
POTENTIAL 

Our work is aimed at obtaining a two-body density- 
dependent potential v p 12 that mimics the three-body po- 
tential. Hence, our starting point is the requirement that 
the expectation values of V123 and of v±2(p) be the same: 



00 

A 



A 



implying in turn (compare to Eqs. ([42]) and ([48f 

E 3 / d ^Vl23 9r23 = E V 12(P) 9l2 ■ 



(50) 



(51) 



A diagrammatic representation of the above equation, 
which should be regarded as the definition of the V12 (p) , 



is shown in Fig. [5] The graph on the left-hand side 
represents the three-body potential times the three-body 
correlation function, integrated over the coordinates of 
particle 3. Correlation and exchange lines are schemati- 
cally depicted with a line having a bubble in the middle, 
while the thick solid lines represent the three-body po- 
tential. The diagram in the right-hand side represents 
the density-dependent two-body potential, dressed with 
the two-body distribution function. Obviously, v p 2 has to 
include not only the three-body potential, but also the 
effects of correlation and exchange lines. 

The left-hand side of Eq. (|5"T]) has been evaluated in 
PI within the FHNC/SOC scheme. Here we discuss 
the derivation of the explicit expression of the two-body 
density-dependent potential appearing in the right-hand 
side of the equation. The procedure consists of three 
different step, each corresponding to a different dressing 
of the diagrams involved in the calculation 

For each of these steps the final result is a density- 
dependent two-body potential of the form 



vi2{p) =^2v p (p,r 12 )O p 2 , 



(52) 



where, depending on the step, the v p {p, ri 2 ) 



— „p 



(p) can 

be expressed in terms of the functions appearing in the 
definition of the UIX potential, the correlation functions 
and of the Slater functions. 



Step I. Bare approximation 

As a first step in the derivation of the density- 
dependent potential one integrates the three-body po- 
tential over the coordinate of the third particle 



(53) 
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Diagrammatically the above equation implies that nei- 
ther interaction nor exchange lines linking particle 3 with 
particles 1 and 2 are included. Only the two-body distri- 
bution function is taken into account in the calculation 
of the expectation value of V123 



(V0 

A 



El 

3! 



E / 



E 

p 



dx 3 V£ 3 



P 

9l2 



(54) 



Note that only the scalar repulsive term and one per- 
mutation of the anticommutator term of the three-body 
potential provide non vanishing contributions, once the 
trace in the spin-isospin space of the third particle is 
performed. 

As shown in Fig QjJJ the contribution of the density- 
dependent potential to the energy per particle of SNM 
and PNM (vip p )/A IS more repulsive than the one ob- 
tained from the genuine three-body potential UIX. Thus, 
the scalar repulsive term is dominant when the three- 
body potential is integrated over particle 3. 



A. Step II. Inclusion of statistical correlations 

As a second step we have considered the exchange lines 
that are present both in 5123 and g\2- Their treatment is 
somewhat complex, and needs to be analyzed in detail. 

Consider, for example, the diagram associated with 
the exchange loop involving particles 1, 2 and 3, de- 
picted in Fig. [6J Its inclusion in the calculation of the 
density-dependent two-body potential would lead to dou- 
ble counting of exchange lines connecting particles 1 and 
2, due to the presence of the exchange operator Pi 2 in 
<7i2. This problem can be circumvented by noting that 
the antisymmctrization operator acting on particles 1, 2 
and 3 can be written in the form 



1 



P12 
(1- 



-Pl3 
PIS" 



" p23 + P12P13 + P13P23 
P23) x (1 - P12) , 



(55) 



in which the exchange operators contributing to the 
density-dependent potential only appear in the first term 
of the right-hand side. 

On the other hand, the second term in the right-hand 
side of Eq. ([55jl only involves the exchange operators 
P12, whose contribution is included in g\i and must not 
be taken into account in the calculation of i>i2(p)- 

Two features of the above procedure need to be clari- 
fied. First, it has to be pointed out that it is exact only 
within the SOC approximation that allows one to avoid 
the calculation of commutators between the exchange op- 
erators Pi 3 and P23 and the correlation operators act- 
ing on particles 1 and 2. The second issue is related to 
the treatment of the radial part of the exchange opera- 
tors. Although it is certainly true that one can isolate the 
trace over the spin-isospin degrees of freedom of particle 
3, arising from P 13 and P23, extracting the Slater func- 
tions from these operators is only possible in the absence 
of functions depending on the position of particle 3 [39j . 




Figure 6. Three particle exchange loop. 



However, this restriction does not apply to the case 
under consideration, as both the potential and the corre- 
lations depend on r*i3 and F23 . As a consequence, retain- 
ing only the Pi 3 and P23 exchange operators involves an 
approximation in the treatment of the the Slater func- 
tions, whose validity has been tested by carrying out a 
numerical calculation. 

By singling out the radial dependence of the exchange 
operators, and by computing the inverse of the operator 
(1 — P12), where pj denotes the spin-isospin part of Pij, 
it is possible to find a "Slater Exact" density-dependent 
potential vf 2 E '{p) whose calculation does not involve any 
approximations concerning the Slater functions. It can 
be rewritten in the form 



S.E. 

J 12 



(P) = 



P 



dx 3 V x 



.,{ 



1 



1 



'12 



P23(^12^13^23 — ^2 



^12^13) + -P13P23 (^12^13^23 



-Pl3(^l2^13'23 — ^13) 
Pl2P23^12^13^23 — 

i?a&)]}, (56) 



Note that in the above equation we have omitted all 
correlations functions, whose presence is irrelevant to the 
purpose of our discussion. The density-dependent poten- 
tial obtained from Eg. ((56)) must be compared to the one 
resulting from the approximation discussed above, which 
(again neglecting correlations) leads to the expression 



,S.A. 



vxr-(p) = I / dx s v 123 (i - p 13 q 3 - p 23 q 3 ) , (57) 

where "S. A." stands for Slater Approximation. We have 
computed {vf 2 E '{p)) and (vf 2 A (p)) for SNM within the 
FHNC/SOC scheme, for both the scalar and the anti- 
commutator terms of the UIX potential. 

The results, plotted in Fig. [71 clearly show that 
Eq. (|57p provides an excellent approximation to the exact 
result for the exchanges of Eq. (|56| . Hence it has been 
possible to use Eq. (f57|) also to compute the contribu- 
tion coming from the commutator of the UIX potential, 
avoiding the difficulties that would have arisen from an 
exact calculation of the exchanges. 

The second step in the construction of the density- 
dependent potential is then 
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Scalar contribution 




0.1 0.2 0.3 0.4 0.5 

P (fm- 3 ) 



Anticommutator contribution 




0.1 0.2 0.3 0.4 0.5 

P (fm" 3 ) 



Figure 7. Contributions of the density-dependent potential to 
the energy per particle (see Eqs. (|56|) and (|57l) ~). arising from 
the scalar term of UIX (upper panel) and from the anticom- 
mutator term (lower panel). 
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Figure 8. Diagrams contributing to the density-dependent 
potential. The dashed lines with diamonds represent the first 
order approximation to g\,J^° i r ij), discussed in the text. Only 
diagrams with at most one operator attached to a given point 
are taken into account. 



v{l{p) = v?i A ip) (58) 

which is a generalization of the bare potential of Eq. (|53| . 

Figure [10] shows that taking exchanges into account 
slightly improves the approximation of the density- 
dependent potential. However the differences remain 
large because correlations have not been taken into ac- 
count. 



B. Step III. Inclusion of dynamical correlations 

The third step in the construction of the density- 
dependent potential amounts to bringing correlations 
into the game. We have found that the most relevant 
diagrams are those of Fig. [5J 



Note that, in order to simplify the pictures, all in- 
teraction lines are omitted. However, it is understood 
that the three-body potential is acting on particles 1, 2 
and 3. Correlation and exchange lines involving these 
particles are depicted as if they were passive interaction 
lines. Moreover, in order to include higher order cluster 
terms, we have replaced the scalar correlation line f?j 
with the Next to Leading Order (NLO) approximation 
to the bosonic two-body correlation function: 

fij 2 -> Qboseinj) = ftj 2 [} + pf df 3 h 13 h 23 ) . (59) 

The full bosonic gboseijij) or gdd(rij) might be used in- 
stead of the NLO approximation. However, including 
higher order terms would have broken our cluster expan- 
sion. The correction to /, c - 2 of Eq. (|59[) . whose diagram- 
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anticommutators of spin-isospin operators, as well as the 

O O = 1 -|- O -O use of suitable angular functions needed to carry out the 

1 2 1 2 integration over r 3 . 
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Figure 9. NLO approximation to the bosonic two-body cor- 
relation function. 



matic representation is displayed in Fig. [9j can indeed 
be considered to be of the same order as the operatorial 
correlations. 

Figure [5] shows that the vertices corresponding to par- 
ticles 1 and 2 are not connected by either correlation 
or exchange lines. All connections allowed by the dia- 
grammatic rules are taken into account multiplying the 
density-dependent potential by the two-body distribution 
function, according to the definition of Eq. ([5Tjl . 

We have already discussed the exchange lines issue, 
coming to the conclusion that only the exchanges Pi 3 and 
P23 have to be taken into account. This is represented 
by the second diagram, where the factor 2 is due to the 
symmetry of the three-body potential, that takes into 
account both Pi 3 and P23. 

The explicit expression of 1*12 (p) obtained including 
the diagrams of Fig. [5] can be cast in the form 



12:! 



9bose l r 13)6W (723) 



NLO 1 



x (1 - 2Pi 3 ^ 3 ) + 4g^ e °(r 13 )f c (r 23 )f(r 23 ) 



a 



■a 



o 

Oh 



> 



a 
o 
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c 
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bo 

a 



o 

Oh 





(60) 



Figure 10. Contributions of the density-dependent potential 
to the energy per particle of SNM (upper panel) and PNM 
(lower panel), compared to the expectation value of the gen- 
uine three-body potential UIX: (Vi23)/A 



where f(r2 3 ) denotes the sum of non central correlations 



/>2 3 ) = E/ p ^ 3 )of J 



(61) 



Note that, in principle, an additional term involving the 
anticommutator between the potential and the correla- 
tion function should appear in the second line of the 
above equation. However, due to the structure of the 
potential it turns out that 



dx 3 {Vl23, /0"23)} 



dx 3 Vi2 3 f(r 2 3) 



(62) 



The calculation of the right-hand side of of Eq. (f6"0"l) 
requires the evaluation of the traces of commutators and 



As for the previous steps, we have computed the con- 
tribution of the density-dependent potential v^ 11 ^ (p) to 
the energy per particle. The results of Fig. [10] demon- 
strate that the density-dependent potential including cor- 
relations is able to reproduce the results obtained using 
genuine three-body UIX to remarkable accuracy. 

To simplify the notation, at this point it is convenient 
to identify 



V 12 (p) = Wia"^) 



(63) 



Note that the above potential exhibits important dif- 
ferences when acting in PNM and in SNM. For example, 
in SNM v p (p, 7-12) 7^ for p = 1, CT12T12, i5i2Ti2, while in 
PNM vP(p,r 12 ) ^ for p = 1, a 12 , S 12 . 
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V. NUMERICAL CALCULATIONS 

A. Variational approach in FHNC/SOC 
approximation 

An upperbound to the binding energy per particle, 
Ey /A, can be obtained by using the variational method, 
which amounts to minimizing the energy expectation 
value (H)/A with respect to the variational parameters 
included in the model. Its cluster expansion is given by 



(H' 

- L -p =T F + (AE) 2 + higher order terms 



(64) 



where Tp is the energy of the non interacting Fermi gas 
and (AE) 2 denotes the contribution of two-nucleon clus- 
ters, described by the diagram of Fig 1111 Neglecting 
higher order cluster contributions, the functional mini- 
mization of (H)/A leads to a set of six Eulcr-Lagrangc 
equations, to be solved with proper constraints that force 
f c and /(p >x ) to "heal" at one and zero, respectively. 
That is most efficiently achieved through the boundary 
conditions [ll], Hlf 



fP(r > dP) = 8 pl , 
dfP(r) 



dr 



\<n< 



= 0. 



(65) 



Numerical calculations are generally carried out using 
only two independent "healing distances" : d c = d p=1 " A 
and di = d 5 ' 6 . 



The energy expectation value (H) /A, calculated in full 
FHNC/SOC approximation is minimized with respect to 
variations of d c , dt, (3 P , and a p . 

To determine the best values of the variational parame- 
ters we have used a version of the "Simulated annealing" 
algorithm |4£j. In metallurgy the annealing procedure 
consists in heating and then slowly cooling a metal, to 
decrease the defects of its structure. During the heating 
the atoms gain kinetic energy and move away from their 
initial equilibrium positions, passing through states of 
higher energy. Afterwards, when the metal slowly cools, 
it is possible that the atoms freeze in a different configu- 
ration with respect to the initial one, corresponding to a 
lower value of the energy. 

In minimization problems the analog of the position of 
the atoms are the values of the parameters to be opti- 
mized, in our case d c , dt, fi p and a p , while the energy 
of the system correspond to the function that has to be 
minimized, that in our case is the variational energy 



Ey — Ey(d c ,d t , f3 p ,a p ) 



(68) 



In the simulated annealing procedure, the parameters 
d c , dt, (3 p , a p are drawn from the Boltzmann distribution, 
cxp(— Ey/T), where T is just a parameter of the simu- 
lated annealing algorithm, having no physical meaning. 

We have used a Metropolis algorithm, with acceptance 
probability of passing from the state s = {d c ,dt, (3 p ,a p } 
to the proposed state s' = {d' c , d' t ,f3' p , a' p } given by 
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Figure 11. Diagrammatic representation of the two-body clus- 
ter contribution (AE)2 of Eq. (|64[) . The thick lines repre- 
sents both the potential and a kinetic contribution, involving 
derivatives acting only on the correlation functions. The ef- 
fect of the other derivatives is included in 7>. 



Additional and important variational parameters are 
the quenching factors a p whose introduction simulates 
modifications of the two-body potentials entering in the 
Euler-Lagrange differential equations arising from the 
screening induced by the presence of the nuclear medium 



5> p ^(r 8J )0£ 
p=i 



(66) 



The full potential is, of course, used in the energy ex- 
pectation value. In addition, the resulting correlation 
functions f p are often rescaled according to 



F, 



p=l 



(67) 



exp 



E(s')-E(s) 



T 



(69) 



By looking at the distribution of the parameters result- 
ing from the Metropolis random walk, it is possible to 
find the values d c , dt, (3 p and a p corresponding to the 
minimum of Ey, e.g. to the maximum of the Boltzmann 
distribution. As the fictitious temperature T is lowered, 
the system approaches the equilibrium and the values of 
the parameters get closer and closer to d c , d t , /3 p ,a p . 

The best solution found during the execution of the 
Metropolis algorithm has been kept. The discrete values 
of the temperature, as well as the numbers of Monte 
Carlo steps for each Tj have been chosen in such a way 
that different executions of the simulated annealing pro- 
cedure produce the same value for d c , d t , fi p and a p . 

A constrained optimization has been performed, by 
imposing the sum rules for the kinetic energy and for 
the scalar two-body distribution function. In particu- 
lar the difference between the Pandharipande-Bethc (PB) 
and the Jackson-Feenberg (JF) kinetic energies has been 
forced to be less than 10 % of the Fermi Energy Tp of 
Eq. (|64| . while the sum rule ([45]) for g c {r 12) has been 
satisfied with a precision of 3 %. 

In our calculations we have optimized the variational 
paremeters for four different Hamiltonians, each corre- 
sponding to different potential terms: Argonne v' s , Ar- 
gonnc i)g + UIX, Argonne v' 6 , and Argonne v' 6 +UIX. 
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Figure 12. Energy per particle for PNM, obtained using the 
density-dependent potential of Eq. (|53|) added to the Argonnc 
v'g (upper panel) and to Argonne Vg (lower panel) potentials. 
The energies are compared to those obtained from the genuine 
three-body potential and from the two-body potentials alone. 



The energy per particle of SNM and PNM computed 
adding to the two-body potentials Argonne v' 8 and Ar- 



gonne v' 6 the density-dependent potential of Eq. (|52j) . 
have been compared to the results obtained using the 
hamiltonian of Eq. ([1]) with the same two-body potentials 
and the Urbana IX three-body potential. The energy as- 
sociated with the density-dependent potential has been 
computed with the same variational parameters result- 
ing from the genuine three-body potential, i. e. no opti- 
mization procedure has been performed for the density- 
dependent potentials. 

Both calculations have been consistently carried out 
within the FHNC/SOC scheme. 

It is worth noting that our simulated annealing con- 
strained optimization allows us to: i) reduce the viola- 
tion of the variational principle due to the FHNC/SOC 



approximation; ii) perform an accurate scan of the pa- 
rameter space. As a consequence, our FHNC/SOC cal- 
culations provide very close results to those obtained via 
Monte Carlo calculations, as shown in Figs. [T2] and [T21 
to be compared with those of Ref. [1] where the agree- 
ment between FHNC and Monte Carlo methods were not 
nearly as good. 



B. Auxiliary Field Diffusion Monte Carlo 
(AFDMC) approach 

In order to check the validity of our variational 
FHNC/SOC calculations, we carried out AFDMC sim- 
ulations dJ for both PNM and SNM. 

The AFDMC method has proved to be a powerful 
approach to deal with large nuclear systems, such as 
medium-heavy nuclei and infinite matter. Using a fixed- 
phase like approximation, AFDMC also yields results in 
very good agreement with those obtained from Green 
Function Monte Carlo (GFMC) calculations for light nu- 
clei 01. 

We have computed the equation of state of PNM and 
SNM using the AFDMC method with the fixed-phase 
like approximation. We simulated PNM with A — 66 
and SNM with A = 28 nucleons in a periodic box, as 
described in [42| and 43 1. The finite-size errors in PNM 
simulations have been investigated in [43j by comparing 
the twist averaged boundary conditions with the periodic 
box condition. It is remarkable that the energies of 66 
neutrons computed using either twist averaging or peri- 
odic boundary conditions turn out to be almost the same. 
This essentially follows from the fact that the kinetic en- 
ergy of 66 fermions approaches the thermodynamic limit 
very well. The finite-size corrections due to the inter- 
action are correctly estimated by including the contribu- 
tions given by neighboring cells to the simulation box[45|. 
From the above results for PNM and those reported in 
[J] for SNM, we can estimate that the finite size errors 
in the present AFDMC calculations do not exceed 1% of 
the energy. 

The statistical errors, on the other hand, are very small 
and in the Figures are always hidden by the squares, 
the triangles and the circles representing the AFDMC 
energies. 



C. PNM equations of state 



In the PNM case (sec Fig. [T^J, the EoS obtained 
with the three-body potential UIX and using the density- 
dependent two-body potential are very close to each 
other. For comparison, in Fig. [12] we also report the re- 
sults of calculations carried out including the two-body 
potential only. In our approximation, with the exception 
of the line with diamonds of Fig. [5J we have neglected 
the cluster contributions proportional to p 2 . One could 
then have guessed that the curves corresponding to the 
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Table I. Values for the saturation densities, the binding energy 
per particle, and the compressibility of SNM obtained from 
the variational FHNC/SOC EoS of Fig. [13] 





Vg + Ul23 


«6 +v(p) 


V' 8 + Vl23 


v'a + v(p) 


Po (fm- 3 ) 


0.17 


0.16 


0.16 


0.15 


E (MeV) 


-11.3 


-11.2 


-10.3 


-10.3 


K (MeV) 


205 


192 


189 


198 



UIX and density-dependent potential would have slightly 
moved away from each other at higher densities because, 
as the density increases, the contributions of higher or- 
der diagrams become more important. Probably, in this 
case a compensation among these second and higher or- 
der terms takes place. 

The density-dependent potential obtained in the 
FHNC/SOC framework has been also employed in 
AFDMC calculations. As can be plainly seen in Fig. 
[T2l the triangles representing the results of this calcula- 
tion are very close, when not superimposed, to the circles 
corresponding to the UIX three-body potential AFDMC 
results. 



D. SNM equation of state 

In the EoS of symmetric nuclear matter, the above 
compensation does not appear to occur, as can be seen 
in Fig. [T3j At densities lower than p = 0.32 fm - , the 
curves resulting from UIX and the density-dependent po- 
tential are very close to one other, while for p > 0.32fm -3 
a gap between them begins to develop. 

The gap is smaller when the two-body potential Ar- 
gonne v' s is used, but the reason for this is not completely 
clear. 

We have computed the saturation density p$, the 
binding energy per particle E(po) and the compressibil- 
ity K = 9p (dE(p)/dp) 2 for all the EoS of Fig. Ql 
The variational FHNC/SOC results are listed in Ta- 
ble HI while those coming from the AFDMC calcula- 
tion with v'q + vi 2 (p) potential are: po = 0.17fm~ , 
E = -10.9 MeV and K= 201 MeV. 

The saturation densities are quite close to the empirical 
value po = 0.16 fm -3 (MeV). For the genuine three-body 
potential this is not surprising, since the parameter Uq 
is chosen to fit the saturation density, as discussed in 
Section [III On the other hand, the fact that the density- 
dependent potential also reproduces this value is remark- 
able and needs to be emphasized. 

The binding energies obtained with Uia(p) are very 
close to those coming from UIX potential, but they are 
larger than the empirical value Eq = — 16McV. 

As for the compressibility, the experimental value K « 
240 MeV suffers of sizable uncertainties. However, also in 
this case the result obtained with the density-dependent 
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Figure 13. Same as in Fig. [12] but for SNM 



potential differs from that obtained with the UIX poten- 
tial by less than 5%. 



VI. CONCLUSIONS 

We have developed a novel approach, allowing one to 
obtain an effective density-dependent NN potential tak- 
ing into account the effects of three- and many- nucleon 
interactions. 

The resulting effective potential can be easily used in 
calculations of nuclear properties within many-body ap- 
proaches based on phenomenological hamiltonians, in- 
cluding the effects of strong NN correlations, which can 
not be treated in standard perturbation theory in the 
Fermi gas basis. Moreover, the derivation of the density- 
dependent NN potential is fully consistent with the treat- 
ment of correlations underlying the FHNC and AFDMC 
approaches. 

While the reduction of n-body potentials to a two- 
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body density-dependent potential is reminiscent of the 
approach of Refs. (Til fl2l ]. our scheme significantly im- 
proves upon the TNI model, in that i) it is based on a 
microscopic model of the three nucleon interaction pro- 
viding a quantitative description of the properties of few 
nucleon systems and ii) allows for a consistent inclusion 
of higher order terms in the density expansion, associated 
with four- and more-nucleon interactions. 

As shown in Figs. [T3] and [T31 the results of calcu- 
lations of the PNM and SNM equation of state carried 
out using the density-dependent potential turn out to be 
very close to those obtained with the UIX three-body 
potential. In this context, a critical role is played by 
the treatment of both dynamical and statistical corre- 
lations, whose inclusion brings the expectation value of 
the effective potential into agreement with that of the 
UIX potential (sec Fig. [TO]) . This is a distinctive fea- 
ture of our approach, as compared to different reduction 
schemes based on effective interactions, suitable for use 
in standard perturbation theory [46|, [47| . 

Using the density-dependent potential we have been 
able to carry out, for the first time, a AFDMC calcula- 
tion of the equation of state of SNM consistently includ- 
ing the effects of three nucleon forces. The results of this 
calculation show that the Ug+UIX hamiltonian, or equiv- 
alcntly the one including the effective potential, fails to 
reproduce the empirical data. 

The FHNC results obtained using the v' 8 potential indi- 



cate that the 5-6 MeV undcrbinding at equilibrium den- 
sity can not be accounted for replacing the v' e with a 
more refined model, such as vi S . Hence, the discrepancy 
has to be ascribed cither to deficiencies of the UIX model 
or to the effect of interactions involving more than three 
nucleons. 

The immediate follow up of our work is the AFDMC 
calculation of the SNM equation of state with the v' a 
potential and the density-dependent potential, which is 
currently being carried out. Further development will in- 
clude a study of the dependence on the specific model of 
three-nucleon force, as well as the inclusion of of four- and 
many-nucleon interactions, whose effects are expected to 
be critical for the determination of the properties of high 
density neutron star matter. 

As a final remark, the effective potential discussed in 
our paper could be easily employed in many-body ap- 
proaches other than those based on the CBF formalism or 
quantum Monte Carlo simulations, such as the G-matrix 
and self-consistent Green function theories [48l - {50j . 
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